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Abstract. The numerical integration of the Nose- Hoover dynamics gives a 
deterministic method that is used to sample the canonical Gibbs measure. The 
Nose-Hoover dynamics extends the physical Hamiltonian dynamics by the addition 
of a "thermostat" variable, that is coupled nonlinearly with the physical variables. 
The accuracy of the method depends on the dynamics being ergodic. Numerical 
experiments have been published earlier that are consistent with non-ergodicity of the 
dynamics for some model problems. The authors recently proved the non-ergodicity 
of the Nose-Hoover dynamics for the one-dimensional harmonic oscillator. 

In this paper, this result is extended to non-harmonic one-dimensional systems. 
It is also shown for some multidimensional systems that the averaged dynamics for 
the limit of infinite thermostat "mass" have many invariants, thus giving theoretical 
support for either non-ergodicity or slow ergodization. Numerical experiments for a 
two-dimensional central force problem and the one-dimensional pendulum problem 
give evidence for non-ergodicity. 
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1. Introduction 

The computation of equilibrium statistical properties of molecular systems is of great 
importance in materials science, computational physics, chemistry, and biology [6l[ll]. 
These equilibrium statistical properties are given by phase space integrals of the form 

{A) = J Aiq,p)dfxiq,p), (1) 

where q = (gi, . . . , qat) G and p = (pi, . . . ,pn) G denote a set of positions 
and momenta and A{q,p) is an observable, a function defined over the phase space and 
related to the macroscopic quantity under study. The computation of integrals such 
as ([1]) is often a challenging problem, especially when the number of degrees of freedom 
is large. 

For molecular systems at fixed temperature 6, the measure dfi is the Gibbs measure 
for the canonical ensemble [6[ [H] 



dfi{q,p) 



exp {-(3H{q,p)) 
exp {—l3H{q,p)) dqdp 



dq dp, (2) 



where H{q,p) is the Hamiltonian of the system and f3 is related to the temperature 6 by 
P = l/{kB9) with ks denoting the Boltzmann constant. We will consider Hamiltonians 
of the general form 

= + (3) 

where M{q) G M^^^ for q G is the generalized mass matrix and V{q) is the potential 
energy. We assume that the generalized mass matrix M{q) G M^^^ is symmetric 
and positive definite, so its inverse M~^{q) G M^^^ exists for all q G and is also 
symmetric and positive definite. 

Many methods have been proposed and utilized to approximate the phase space 
integral ([1]), including methods based on stochastic or deterministic dynamics for {q,p). 
If the dynamics is ergodic with respect to the measure dfi given by ([2]), then the phase- 
space average ([1]) is equal to the time average 

[ Aiq,p)dfxiq,p)= hm ^ [ A {q{t) , p{t)) dt (4) 

over a trajectory {q(t),p{t))t>o- Thus, the time average can be approximated by 

lim -/ A{q{t),p{t))dt^ Jim — J^MQe^Pi), 

where {qe,Pe)e>i is a numerical solution of the chosen dynamics. 

In this paper, we investigate the deterministic dynamics known as Nose-Hoover 
dynamics [7] , which is still widely used although variants have been developed with the 
goal to improve its efficiency and overcome its deficiencies [T2| [T6| [3l [11] . This dynamics 
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has been first proposed in the form of a Hamiltonian dynamics on an extended phase 
space [15], the Hamiltonian being chosen such that the marginal distribution of its 
microcanonical density is the canonical Gibbs density for the physical variables. The 
Nose-Hoover dynamics is then constructed by rescaling time and momentum to obtain 
a non-Hamiltonian dynamics with physical time and momentum [7]. 

Stochastic dynamics (such as the Langevin equation, or the recently proposed 
Hoover-Langevin method [10]) can also be considered. See [1] for a review of sampling 
methods of the canonical ensemble, along with a theoretical and numerical comparison 
of their performances for molecular dynamics. 

The equality (jlj) relies on an ergodicity condition. This condition has been 
rigorously proven neither for the Nose-Hoover dynamics, nor for any other deterministic 
method commonly used in practice. In fact, there is numerical evidence that shows that 
the Nose- Hoover method is not ergodic for some systems [TJ [I2l [16], including the one- 
dimensional harmonic oscillator. In [9], we have rigorously analyzed the dynamics in 
this special case, and indeed proven the non-ergodicity, for some regime of parameters. 

In this article, we study more general systems. After briefly recalling the Nose- 
Hoover equations (see Section [2]), we first consider a class of multidimensional systems 
(see Section [3]). Taking the limit of an infinite thermostat "mass" in the Nose-Hoover 
equations, we formally obtain an averaged dynamics, for which we prove the existence of 
many invariants. These theoretical results are illustrated by numerical simulations of a 
specific system (see Section H]). We numerically observe that, for finite thermostat mass, 
these invariants are of course not exactly preserved, but still remain close to their initial 
value. This prevents the Nose- Hoover system from thermalizing. In Section [5l we next 
turn to the one-dimensional case, for which we obtain stronger results. We first prove 
non-ergodicity of the Nose- Hoover dynamics, when the mass of the thermostat is large 
enough (see Section [5?T1) . Our method extends the one we used to study the harmonic 
oscillator case [9] . Section 15.21 describes an example of such a one degree of freedom 
problem. Again, numerical simulations illustrate the obtained theoretical results. 



2. Nose-Hoover dynamics 

The Nose-Hoover dynamics involves the physical variables q and p and one additional 
scalar variable, ^, which represents the momentum of a thermal bath exchanging energy 
with the system. The differential equations are: 

P = -^-QP=-^nq) ^ qP, (5) 

T 1.x A^ 

e = P^M~\q)p~-, 

where 1 denotes the time-derivative. The parameter Q represents the mass of the 
thermostat; it is a free parameter that the user has to choose. 
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We recall that invariant measures p{z) dz for a general dynamical system 
z = f{z) 

are determined by the equilibrium equation 
div(p(^)/(^)) = 0. 

It can be verified by direct computation that the dynamics ([5]) preserves the measure 



dfiNH = exp 



dq dp d^ (6) 



by using the fact that the kinetic energy is quadratic in p. 

If the dynamics ([5]) is ergodic with respect to (i/^NH, then, by integrating out C,, we 
have that the dynamics {q(t),p(t)) is ergodic with respect to the Gibbs measure. In 
this case, the time-average of a function A{q,p) along a typical Nose- Hoover trajectory 
provides an estimate for the space-average of A with respect to Gibbs measure. 
Unfortunately, the system is generally not ergodic. In [9j we proved non-ergodicity 
in the case of the one-dimensional harmonic oscillator. Our aim here is to study more 
general systems. 

3. Systems with first integrals 

In this section, we show how the presence of additional integrals for a Hamiltonian 
system can impede ergodization of the Nose-Hoover dynamics. 

3.1. Homogeneous integrals 

Consider a Hamiltonian system 

. _dl£ . _ _dH_ 
^ dp ' dq ^ 

for energy ([3]) which admits a first integral other than H itself. This means that there 
is a smooth function F{q,p) whose Poisson bracket with H vanishes, i.e., 

{H, F} = HjF, - HjF, = 0. (8) 

If F is a homogeneous function of the momentum variables, then it gives rise to a first 
integral of the Nose-Hoover system. 

Theorem 3.1 If F{q,p) is a first integral of ^ which is homogeneous of degree k with 
respect to the momentum variables, p, then 

G{q, P^O = ^ + H{q,p)~jj: In \F{q, p) \ (9) 
is a first integral of the corresponding Nose-Hoover system 

The proof is a simple computation using ([8]), ([5]) and the fact that F^ p = kF. 

Of course, the existence of such an integral immediately gives non-ergodicity of the 
Nose-Hoover system with respect to ([6]). For a simple example of a system admitting 
such a homogeneous integral, see Section HI 
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3.2. Completely integrable systems and action-angle variables 

We now assume that the Hamiltonian dynamical system (I7j) is completely integrable, 
i.e., the system admits N independent first integrals which commute in the sense that 
the Poisson brackets of any two of them vanish [1]. The rest of this section is devoted 
to showing that these integrals, even if they are not homogeneous, have a deleterious 
effect on the ergodization of the corresponding Nose-Hoover system. 

The non-degenerate level sets of the integrals are A^-dimensional manifolds and 
if they are compact then their connected components are diffeomorphic to the N- 
dimensional torus T^. Moreover, such a torus has a neig hborhood [/ C x in 
which one can introduce symplectic action-angle variables. More precisely, there exist 
angle variables 6 G T^, action variables a E D G M^, and a symplectic diffeomorphism 
ip : U ^ X D which transforms ([7]) to the form 

e = u{a), d = 0. (10) 

Here D is an open subset of M^. Equivalently, the action-angle Hamiltonian H{6, a) = 

dHiO ci] 

H(ip~^(9, a)) is independent of 6 and — ^ — = uj(a). 

oa 

In what follows, it will be convenient to define the angle mapping ipi{q, p) G M.'^ 
and the action mapping ip2{(l, V) ^ by ipi^q^ p) = {ipi{q, p), 4'2{q, v)) ^ K^^- We will 

also use the abbreviated notation = ^ G M^""^ and ^2^, = ^ G R^^^. We 

aq op 

then denote the Jacobian of ip by 



diipi d2ipi 
diip2 d2ip2 



We denote the inverse mapping of ip{q.,p) by (j){6,a). The matrix (Dip) = (D0)" 
has a simple form since is symplectic: 



where 




j = 

The diffeomorphism transforms ([7]) to ffTOj) by the chain rule 

= ( ! I • (12) 





Since ip is symplectic, the dynamics fllUj) is obtained from the Hamiltonian H{6, a). 
Hence, 
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3.3. Recasting the Nose-Hoover dynamics 



We now multiply the Nose- Hoover equations ([5]) for q and p by D^^ to obtain from (fT2l) 
and ffT3l) that 



As we are interested in the regime Q ^ 1, we rescale by 

Using the symplectic property f|TT]) of ^z^, we see that the Nose-Hoover equation 
be then given in the scaled angle-action variables by 



1'^ 




( ^(«) ^ 


/ 


a 







+ e 






^ ) 





a(c?20i)^02 
-a(9i0i)'^02 



()i^M-i(0i)02 



can 



(14) 



3.4. Averaging the fast variables 



We next apply the averaging method to obtain an approximate system which does 
not involve the fast variables. Rigorous results about averaging for Hamiltonian 
systems with several degrees of freedom are fraught with technical difficulties (see for 
example [2]). These arise from the fact that uj{a), the frequency vector of the fast 
angles, experiences resonances of the form k ■ uj{a) = 0, k & Z^, for certain values of 
the action vector a. In fact, these resonant actions are generally dense in the action 
domain D. In spite of this, the averaged differential equations often provide a useful 
first approximation to the behavior of the slow variables when e is small. 
In our problem, the averaged system for the slow variables is given by 

a = —a S{a), 
a = k{a), 

where 



(15) 



Sia) = ((9i0i)^02)(«)= / {d,<P^fi9,a)M0,a)de, 
k{a) = (0^M-i(0i)02)(a)-|^ 

0^^, a)M-i(0i(^, a))02(^, a) d9 - ^. 



(16) 



We next show that, with an additional assumption on the action-angle mapping (f), 
we have 

S{a) = a. (17) 

Recall that a map {q,p) = (f){6,a) is symplectic if it preserves the canonical differential 
two-forms, i.e., (J)* dpi A dqi) = Y^daiAdOi where 0* indicates the pull-back, meaning 
that we write pj, g,, dpi, dqi in terms of the 6, a variables. It follows from this that the 
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difference of the corresponding canonical one- forms 0* {p^ dq) — a'^ dO is a closed one- 
form on X D. We will call (f) exact symplectic if this closed one-form is exact, i.e., 
if 

(f)*{p^ dq) =a^de + dF{e, a) (18) 

where F{6, a) is a real-valued function on x D. This stronger condition holds for 
the action-angle coordinates associated to many well-known integrable systems. 

As an example, consider the one- dimensional harmonic oscillator H{q,p) = |(p^ + 
q"^) which is a completely integrable system with = 1 degrees of freedom. In any 
annulus of the form < /ii < (p^ + g^)/2 < /i2, we can introduce action-angle variables 
{9, a) such that 

g = v^2acos^, p = —V2asin6. 

For (f>{9,a) = {q,p), we have 

0* (p dq) = — sin ^ cos 9da + 2a sin^ 9d9 



and so 



where 



4>*{pdq) — ad9 = — sin 9 cos9da + a(2 sin^ 9 — l)d9 = dF 



F = —a sin 9 cos 9. 

It turns out that the action-angle variables constructed according to the usual 
method of Arnold [1] always have this exactness property. To see this, recall that 
in Arnold's method, the tori given by fixing the A^-independent integrals of motion 
are parametrized by angle variables 9 = {9i, . . . ,9]\f) derived from the A^ commuting 
Hamiltonian flows defined by the integrals. Then the action variables are given by 

dq (19) 



where 7^ is the curve in the torus defined by holding 9j = const for j ^ i and letting 9i 
run over [0, 1]. The integral depends on which torus is considered, i.e., it is a function 
of the A^ first integrals. The usual proof shows that the map {q,p) = (f){9,a) is defined 
and symplectic on some domain of the form x D. It follows that 

u = (p* (p^ dq) — ad9 

is a closed differential one- form on x D. Showing that is exact symplectic amounts 
to showing that u is exact. For this, it suffices to check that its integral around any 
closed curve vanishes. In fact, since u is closed, it suffices to check the curves of the 
form Ci = Ti X {ao},ao G D where Ti = {9 : 9j = const, j 7^ i}. For such a curve, we 
have 

z/ = / (p* {p dq) — a d9 = / pdq— / a^i d9i = aQi — a^i = . 
c, Jc, Jo 



Here we used the fact that under the curve Ci maps to the curve 7^ used in (jin 
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To prove (IT7|) under the exact symplectic assumption fllSp . we first note that 

(j)*{p^ dq) = (f)2d(f)i = 02 (^^101 d9 + (920i da). 

Hence, ( fTSi) reads 

0^((9i0i rf^ + (9201 rfa) =a^ de + dF{e, a). 

For any j = 1, . . . , A^, we can integrate both sides with respect to 6j, along the circular 
loop Cj as in the last paragraph. Using the periodicity of F in 6j, we obtain 

We can then further integrate fl20l) over the angles 6k for k ^ j to obtain that 



i=l 

forj = l,...,iV. 

5.5. Fzrsi integrals of the averaged Nose-Hoover equations 

A direct calculation shows that a set of independent first integrals for the averaged 
Nose-Hoover equations 



are given by 



a = —aa, 
a = k{a), 



Gi{a, a) = — , 2 = 1, . . . ,iV - 1, 



(21) 



2 k(s^,...,S^,s) (22) 

Gat (a, a) = h / — ^ —ds. 

2 J s 

To prove that G]\f{a, a) is a first integral, it is helpful to use the fact that Gi{a, a) = 
ai/ajq for z = 1, . . . , — 1 is a first integral. 

We summarize the result of this section in the following theorem. 

Theorem 3.2 The averaged equations for the Nose-Hoover dynamics for a completely 
integrable Hamiltonian system has N independent first integrals. 

To the extent that the averaging method applies, we expect that Gi{a{t), a{t)) 
evolves slowly for small e and so the sampling of the Gibbs measure is slow even if the 
dynamics is ergodic. We will verify this numerically in an example in the next section. 
It turns out that Gi{a{t), a{t)) remains quite close to its initial value for fairly large 
values of e as well. 
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4. A central force problem 

We consider here a two degrees of freedom system to illustrate the theoretical results 
obtained in the previous section. We work with the Hamiltonian ([3]) with N = 2, 
the identity mass matrix M{q) = I2, and a potential V^d^'l) which depends only on the 
distance to the origin. The Hamiltonian system ([71) admits two first integrals, the energy 
H and the angular momentum 

L = qiP2 - q2Pi, 

which satisfy {H, L} = 0, and whose gradients are linearly independent, except for values 
of H and L satisfying a condition of the form f{H, L) = for some function /. Hence, 
this system is completely integrable. Assume that V{r) +00 as r ^ +00. Then level 
sets of H are compact, hence the level sets {{q,p) G M"^; H{q,p) = h,L{q,p) = i} are 
also compact, hence there exists action- angle variables for this system. 

To describe the action variables, first introduce polar coordinates (r, 0) in R^. The 
angular momentum is 

L = r20. 

Fixing a value for L, we have a reduced Hamiltonian system for the radial variables 
{r,Pr), where p^ = with Hamiltonian 

1 1 

HL{r,Pr) = -pI + --^ + V{r). 

This reduced system has one degree of freedom and can be understood by the usual 
phase-plane method. Since V{r) +00 as r ^ +00, the level curves 

= {{r.p,) : HLir,Pr) = h}, L 0, 

generically consist of one or more simple closed curves. For the unreduced system, where 
we remember the angle 0, each such curve becomes an invariant torus T(h,L)- It can be 
shown that the action variables assigned to such a torus by Arnold's procedure are as 
follows: ai{h,L) is the area in the {r,pr) plane enclosed by the simple closed curve of 
the reduced system. 



ai{h, L) = Pr dr, 



and a2{h,L) = L, the angular momentum. Note that ai{h,L) is easily computable by 
standard numerical integration schemes. 

Since L is homogeneous of degree A; = 1 in the momentum variables. Theorem 13.11 
gives a first integral for the Nose-Hoover system: 

G{q, P^O = ^ + Hiq,p)-^ In |L(g, p) \ . (23) 

In addition. Theorem 13.21 provides additional integrals for the averaged Nose-Hoover 
equations, in particular the ratio of the action variables 

Qi(g'P) _ (^iiH{q,p),L{q,p)) 

a2{q,p) L{q,p) 
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To the extent that the averaging method apphes for this two-degrees of freedom problem, 
this ratio should evolve only very slowly when Q ^ 1. In the sequel, we present some 
numerical simulations showing first that, when Q ^ 1, this is indeed the case, and 
second that, for Q = 1, such a behaviour persists to some extent. 

Consider the example with potential 

V{r) =r^ + 

with an initial condition (gcPc^o) such that L(go,Po) 7^ 0. We compute the trajectory 
of the Nose- Hoover dynamics ([5]) with the algorithm proposed in [13j. On Figured], we 
plot G{q{t),p{t),^{t)) and Gi{q{t),p{t)), where G and Gi are defined by ([23]) and (IST 
for Q = 100. We indeed observe that G is preserved, whereas Gi evolves slowly. 



1.04 




50000 



t 



Figure 1. Plot of G{q{t),p{t),^{t)) and Gi{q{t),p{t)) (renormalized by their initial 
value) along the trajectory of ([5]), for Q — 100 (/3 = 1, initial condition q = (0;0.5), 
p=(-1.5;1.5),e = 0). 



We now consider the value Q = 1 and plot the same quantities as above on Figure [2] 
Again, G is preserved, whereas Gi evolves in a band which is still quite narrow, even 
for this small value of Q. 

Let us now derive another quantity, which does not behave as well as Gi for large 
Q, but happens to behave in a better way for small (Jil. From the Nose-Hoover dynamics 
([5]), we compute that 



-eaL, 



H 



f 

i rri rrn rri 

j^P P = ~£otp P = —^ot 02 [6, a) (j)2{6, a). 



(25) 
(26) 



I We have a clear understanding of why Gi behaves better than this quantity when Q ^ 1. However, 
the situation for Q = 1 is less clear. 
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10000 20000 30000 40000 50000 

t 

Figure 2. Plot of G{q(t),p{t),^{t)) and Gi{q{t),p{t)) (renormalized by their initial 
value) along the trajectory of ([5]), for Q = 1 (/? = 1, initial condition q = (0;0.5), 
p-(-1.5;1.5),e = 0). 



Since 6 are fast variables whereas a, a and H are slow ones, we again formally use the 
averaging method on fl26l) and consider the dynamics 

H = -ea ko{a) (27) 

with 

ko{a) = I (t)l{e,a)(t)2{e,a)de 



JT2 

(note that (!25l) does not depend on the fast variables 9). Now recall that the action 
variables a are functions of H and L. The equation (l27j) hence reads 

H = -ea ko{H,L). 
The averaged system is thus 
L = —aL, 

H = -ako{H,L), (28) 
d = ko{H,L) -2l3-\ 



oP- 2 

E{E,L,a) = — ^H--\^\L\ (29) 



Note that 



is a first integral of the above system. It is just the analogue of (1231) being a first integral 
for the Nose-Hoover system (see Theorem 13.11) . 

On Figure[3], we plot the function H ^ ko{H, L), for several values of L. We observe 
that ko{H, L) is almost a constant with respect to L, and can hence be approximate(i§| 
by a function k^^{H). 

§ In practice, we have considered several energy values Hi, and for each Hi, we have considered 
several configurations (qi_j,pi_j) with energy Hi and angular momentum Lij. We next have computed 
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_L 
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2 
H 
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Figure 3. Plot oi H ^ ko{L,H). For each value of H, we have considered several 
values of L. 



We hence approximate (128|) by 
L = —aL, 

H = -ae^^^H), (30) 
a = A;^PP(/7) - 2/3-1. 

Now, it is natural to introduce the variable r defined by 



and its reciprocal H{t), such that (!30|) reads 
L = — aL, 

f = — «r, (31) 
a = ko'''{H{T))-2(3-\ 

This system is in the form fl2Tl) . Its two first integrals are 

Ei{L,t) 



T 

L 



and 



E2iT,a) = —+ J —ds 

= ^ / " ^ ^ - -Inr 

2 J s (3 

2 



The first invariant E given by (!29i) is not independent from Ei and E2: E2 = 
E -2(3-^\n\Ei\. 

ko{Hi, Li,j) by averaging p(t)-^p(t) along a constant energy trajectory. Averaging these ko{Hi, Lij), we 
obtain kQ^^{Hi), which next leads to kQ^^{H) for any H by piecewise linear interpolation. 
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We now consider the same trajectories of the Nose-Hoover dynamics that we 
considered on Figures [1] and [2l and we plot 

E^iq,p) = E^ r(if(g,p))) . 

We see on Figure H] that, for Q = 100, this quantity is almost preserved, and that, even 
for Q = 1, it remains close to its initial value. 



1 I I 

Q = l 




10000 20000 30000 40000 50000 

t 

Figure 4. Ei{q{t),p{t)) (renormalized by its initial value) along the trajectory of (O, 
for Q = 1 and Q = 100 (/3 = 1, initial condition q = (0; 0.5), p = (-1.5; 1.5), £, = 0). 



We finally consider an initial condition such that L{qQ, pq) = 0. Along the trajectory 
of (EI), we have L{q(t),p(t)) = by ( 125|) . hence Ei is not defined. On Figure [5l we plot 

E,{q,p,0=E2(^TiHiq,p)),-^ 

along two trajectories, obtained with the same initial condition and the choices Q = 100 
and Q = 1. We again observe that E2 is almost constant for Q = 100, and that it remains 
close to its initial value for Q = 1. 

On FigureEl we plot the energy H {q{t) , p{t)) along the same trajectory (for Q = 1). 
We see that values h < 1 are not sampled. However, there exist {q,p) G such that 
L{q,p) = and H{q,p) is as close to as wanted. Hence, the trajectory only samples 
a strict subset of the level set {{q,p); L{q,p) = 0}. 



5. Systems with one degree of freedom 

Consider a Hamiltonian system of the form ([3]) with = 1 and M{q) = 1. All 
such systems are completely integrable since H itself provides the required integral 
of motion. Suppose there is an interval of energies / = [hi, such that the level curves 
M{h) = {{q,p) : H{q,p) = h},h G /, are all simple closed curves (one-dimensional 
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t 



Figure 5. E2{q{t) , p{t) , ^{t)) (renormalized by its initial value) along the trajectory 
of dSI), for g = 1 and Q = 100 (/3 = 1, initial condition q = (-0.5; 0.5), p = (-1; 1), 

e-0). 




50000 



t 



Figure 6. H{q{t),p{t)) along the trajectory of ([5]), for Q — 1 {P = 1, initial condition 
g= (-0.5;0.5),p= (-1;1), e = 0). 



tori) which are non-degenerate in the sense that the gradient of H does not vanish. 
Then the plane region U = {{q,p) '■ hi < H{p,q) < is diffeomorphic to an annulus, 
and we can introduce action-angle variables (a, 6) in U, and an exact symplectic map 
4>{9,a) = {q,p), as in Section [3l2l 

From Theorem 13.21 we have = 1 first integrals for the averaged Nose-Hoover 
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equations. Let us now rewrite this first integral more explicitly. In view of fl22l) . we have 



«2 



G(a,«) = y + iy(a) (32) 



W{a) = ^ ds, (33) 



with 



where, in view of ( 1161) . k is given by 

k{a) = [ ^1(9, a)de-]-. (34) 

As in the multidimensional case, this integral prevents rapid ergodization, at least 
for small e. But in the one-degree of freedom case we go further and identify conditions 
on H which rigorously imply non-ergodicity. The method is essentially the one used in 
[9] where we treated the harmonic oscillator. Namely, the integral G leads to invariant 
tori of the averaged system which, under certain assumptions, persist for small values 
of e. 



5. 1 . Proof of non-ergodicity 

We will apply a KAM theorem to the Nose-Hoover equations, in the formulation f[T^ . 
Let us introduce the Poincare return map, Pg, of the system (fT^ to the Poincare section 
defined by = mod 1. It is convenient to rescale time by uj{a), so that the return 
time when e = is L This just alters the parametrization of the solutions so that the 
return time to the Poincare section is 1. 

Since there is only one degree of freedom, the averaging method can be rigorously 
justified. Indeed we can eliminate the fast angle 9 of (fT4l) by a change of variables. 
We construct functions g{d,9,a) and h{d,9,a) and corresponding new variables (a, a) 
defined by 

a = d + eg{d,9,a), 
a = a + eh{d,9, a), 

so that in the new variables (d, a), the dynamics IHM is given (after replacing (d, a) by 
(a, a)) by 

9 = u;{a) + 0{e), 

a = -£aS{a) + 0{e^), (36) 
a = ek{a) + 0(£:^), 

where S{a) and k{a) are the averages (fT6ll . 

In view of fl35|) and fl36l) . the Poincare map Pe(d, a) is an 0(£:^) perturbation of the 
time e advance map of the averaged system f|T5|) . for which G defined by fl32|) is a first 
integral. So we now make some assumptions about the level curves of G. 
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Recall that we are working in a region U of the (g,p)-plane defined by an interval 
of energies / which corresponds to an interval of actions J = [01,02]. We assume that 
W{a) has at least one local minimizer ao in J. We have 

ao 

hence k^ao) = and the point P = (ao,0) is an equihbrium point for (fTSl) . The parts 
of the level curves of G which are near P are simple closed curves around P in the 
(a, a)-plane. 



Remark 5.1 If a^ is a local minimizer ofW, then k{ao) = 0, hence j (f)2{d, a^) dO 



T 



Let Go = ^(00, 0) = W{ao) be the value of the integral G at the equilibrium 
point P = (ao,0). Choose constants Gi and G2 such that Gq < Gi < G2 < 
min(G(ai, 0), G(a2, 0)) = mm{W{ai), W{a2)) and let K = [Gi, G2] (see Figure[7]). Then 
the level curves {(a, a); G{a,a) = a}, where cE K, have connected components which 
are simple closed curves near P. The union of these components for c G -ft' forms a 
region D near P which is diffeomorphic to an annulus. 



W(a) 



Gi 
Go 



Figure 7. Schematic representation of tlie interval K = [Gi, G2]. 

The variable G defines a natural action variable in D. We construct the 
corresponding angle variable by following the same method as in [9]. Let Ti{g) 
denote the period of the periodic solutions of (fT5l) which corresponds to the level curve 
G = g & K . The averaged differential equation ( ITSl) becomes 

= l/UG), 

G = 0. ^'^^^ 
In these coordinates, the time e advance map takes the form (0, G) ^ (01, Gi) where 

01 = + e/Ti(G), ^gg^ 
Gi = G . 

Call this map Q,(0, G). Then P,(0, G) = Qe{(j), G) + 0{e^). 

Theorem 5.1 Suppose the period function Ti{G) is not identically constant on the 
interval K . Then, for e sufficiently small, the Poincare map P^ has invariant circles in 
the region D and so the Nose-Hoover system is not ergodic: trajectories {q{t),p{t)) that 
solve ^ are not ergodic with respect to the Gibbs measure 
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Proof: As in [9J, we apply Moser's twist theorem to the Poincare map P^. The details 
are similar to those in |9], so we only sketch the argument here. 

The fact that the Nose-Hoover differential equation preserves the invariant measure 
([6]) implies (as in [9]) that P^ preserves an invariant measure in the (a, a) plane. It 
follows that the maps Pe have the curve intersection property. The hypothesis on Ti{G) 
guarantees that, making K and hence D smaller if necessary, we may assume that either 
T[{G) > or T[{G) < throughout D. This means that there are many invariant circles 
in D for which the rotation number under is Diophantine. Moreover, the required 
twist condition holds. Moser's theorem guarantees that such invariant circles perturb 
to nearby invariant circles for e sufficiently small. 

We now show that the existence of these invariant circles implies non-ergodicity 
with respect to the Gibbs measure. First, note that on a level curve }A = 
{{a, a); G{a,a) = c}, where c G -ft', we have that W{a) is bounded from above, since 
W{a) < G{a,a) = c. In view of the choice of K (see Figure [7]), this implies that 
a is lower and upper bounded. Since a'{h) is positive and bounded away from for 
h E I = [hi,h2], this hence shows that H{q,p) is lower and upper bounded (that is, 
\H{q,p)\ is bounded) on the invariant circle of Qe- As a consequence, \H{q,p)\ 
is bounded on the nearby invariant circles of P^. Hence, the trajectory of ([5]) does 
not sample values of H{q,p) larger than some threshold. This is a contradiction with 
{q(t),p(t)) that solves ([5]) being ergodic with respect to the Gibbs measure ([2]). □ 

Because of the complicated series of coordinate changes leading from the original 
Hamiltonian system to the averaged system, it is not easy to state simple conditions on 
the original potential function V{q) which guarantee that the period function Ti{G) is 
not constant. An equilibrium point surrounded by periodic orbits of constant period is 
called isochronous and various criteria for isochronicity have been given. Our problem 
can be reduced to a Hamiltonian case for which a simple criterion can be stated. 

To carry out the reduction, replace a in (fT5l) (that is, (12T1) ) by a = ln(a/ao), where 
Oq is a local minimizer of W (see Figure [7]). The differential equation fl2Tl) becomes 

a = -a, a = U'{cr), (39) 

where U{a) = W{aoexpa). Except for a reversal of time, this is a classical Hamiltonian 
system with Hamiltonian G{a,a) = o? jl + U{a). It has an equilibrium point at the 
origin (a, a) = (0, 0). 

Now [8] discusses the problem of recovering the potential of such a system from its 
period function (see also [5J). Let Go = ^(0) be the energy level of the equilibrium point 
at the origin. For G > Gq let L{G) be the width of the potential well at energy G, i.e., 
L{G) = o"2(G) — cri(G) where (Ti{G) are the two roots of U{(j) = G near cr = 0. Then 

Ti(G) is constant if and only if L(G) = — a/2(G — Go). This is just the formula for the 

71 

width of the quadratic potential well associated to a harmonic oscillator of period Ti. 
Clearly this is highly exceptional and is easy to rule out, at least numerically. 

Another way to show that Ti{G) is non-constant is to observe that the constancy 
of the period implies that the family of periodic orbits surrounding the equilibrium 



Non-ergodicity of Nose-Hoover dynamics 



18 



point must fill the entire plane If this were not so, then there would be another 
equilibrium point on the boundary of the maximal family which would force Ti{G) —>■ oo. 
For example, for certain values of /?, the pendulum equations (see Section 15.21) lead to 
an averaged system with more than one equilibrium, and this immediately implies that 
Ti{G) is non-constant. 

5.2. The simple pendulum problem 

We consider here the numerical example of a simple pendulum whose potential energy 
is given by 

V{q) = — cos q. 

We reduce q modulo 2n. By construction, the energy satisfies h > —1. The phase 
portrait is shown on Figure [HI The above assumptions are satisfied for energies in the 
interval / = [hi, /i2], with — 1 < /ii < /;,2 < 1, or 1 < /ii < h2. 




Figure 8. Phase portrait of tlie simple pendulum. 



First, we numerically compute a{h) defined by (fT9l) . Note that, in this one- 
dimensional setting. 



a{h) = pdq, 

J M{h) 

where the line integral is taken in the direction of the Hamiltonian flow, and M{h) = 
{{q,p) G : H{q,p) = h}. We also compute 

Jj 

which is independent of P and satisfies k{a) = ko{a) — (3~^, with k defined by ( IMI) . In 
practice, /cq is computed using the fact that, for any energy level h, 



1 f 

A;o(a(/i)) = lim -/ p^{t)dt, 



where {q{t),p{t)) solve the Newton equations of motion for the pendulum at the constant 
energy h. Results are shown on Figure [91 
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-1 -0.5 0.5 1 1.5 2 -1 -0.5 0.5 1 1.5 2 

h h 

Figure 9. Numerically computed values of a{h) and ko{a{h)) (see text). 



We have seen that the action values a such that ko{a) = that is k{a) =0, play 
an important role (see Remark l5. II) . In view of Figure[9], we see that, when j3~^ < (3~^ for 
some threshold /?c, then the equation kQ^a) = [3~^ has three solutions. When (3~^ > 
then the equation fco(a) = has a unique solution. In what follows, we detail the 
numerical results obtained with the choice (3 = 1, which corresponds to the first case. 
Similar results have been obtained for choices of (3 corresponding to the other case. 
Hence, the conclusions that we draw here are by no means restricted to the case (3 = 1. 

The function W{a) defined by fl33|) is shown on Figure [10] for the choice (3 = 1. 
This function has two local minimizers, ai ^ 7.6 and a2 ~ 16.17. Following the above 
theoretical analysis, we work close to one of them. We have chosen to work close to ai. 




a 



Figure 10. Numerically computed values of W{a) for (3=1. 

On Figure [TTl we plot the trajectory of the averaged dynamics f|T5l) for different 
initial conditions. These trajectories have been computed with the Symplectic Euler 
algorithm used on the Hamiltonian formulation (!39|) . As expected, the trajectory is 
a simple closed curve around the equilibrium point (ai,0), that corresponds to a level 
curve of G. These curves are also invariant curves of the map Qg defined in the previous 
section (see map fl38l) ). 
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We now study how these curves persist upon perturbation. We recall that the 
Poincare return map Pg of the Nose-Hoover dynamics f[T^ on the section ^ = mod 1 
is a perturbation of Qs- Results for the Poincare return map of the dynamics f|T^ are 
shown on Figure [T^ for Q = 10^ (that is, e = y/To x 10~^), and on Figure for Q = I 
(that is, e = 1), for the same initial energies as for Figure [TTl These Poincare return 
maps have been computed using the fact that the section = mod 1 corresponds to 
the section g = mod 2tt. We see a good agreement between Figures [TT] and [T^l The 
presence of invariant circles on Figures [12] and [13] shows that the system f [T4|) seems to 
have invariant curves, for Q = 10^ and Q = 1. 

Note that Theorem 15.11 which states the non-ergodicity of the Nose-Hoover 
equations, relies on the important assumption that the period Ti{G) of the averaged 
equations is not constant. This holds true for the pendulum case, in view of the 
discussion at the end of Section 15. 1[ This is also confirmed by numerical computations 
of Ti(G') (see Figure [Tl). 

Let us now look at another criterion for ergodicity, namely what energy values are 
sampled. We see on Figure [I3] that small values of a are not sampled: we have a > 6 
for the three initial conditions that we considered. In view of Figure [9], this corresponds 
to small values of H not being sampled. On Figure [TS] we plot the physical energy 
H{q(t),p{t)) along the trajectory of dHl), for the value Q = 1, and the initial condition 
q = 0, p = 1.5, ^ = 0, that corresponds to the initial value a(0) = 7.72 that we studied 
on Figures [HI [12] and [13] (results are the same for other initial conditions). We see 
that H > —0.4. If the dynamics (HM was sampling the canonical measure, then all 
values of H would be attained. In particular, the smallest values H ^ —1 would be the 
most frequent ones. Indeed, from the Gibbs measure Q, we compute the probability 
distribution function of the energy, which reads p{h) = z^^ exp{—f3h) a'{h), where z 
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is a normalization constant. For the pendulum (h) is close to a constant (see 

Figure E]), hence the smallest values of h are the most frequent ones. Hence, it seems 
that f|T^ is not ergodic with respect to the canonical measure, even for the value Q = 1. 
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Figure 15. Energy H{q{t),p{t)) along the trajectory of ([5]), for Q — 1 and /3 — 1, and 
the initial condition q — 0, p — 1.5, £, — (that is, a(0) — 7.72). 
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